Broad histogram relation for the bond number and its applications 
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We discuss Monte Carlo methods based on the cluster (graph) representation for spin models. 
We derive a rigorous broad histogram relation (BHR) for the bond number; a counterpart for 
the energy was derived by Oliveira previously. A Monte Carlo dynamics based on the number of 
potential moves for the bond number is proposed. We show the efficiency of the BHR for the bond 
number in calculating the density of states and other physical quantities. 
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I. INTRODUCTION 

The development of new algorithms for the Monte 
Carlo simulation is important to overcome the problem of 
slow dynamics. We may classify such attempts into two 
categories. The first category is the extended ensemble 
method; one uses an ensemble different from the ordi- 
nary canonical ensemble with a fixed temperature. The 
multicanonical method BL 0, the simulated tempering 
H , the exchange Monte Carlo method , the broad his- 
togram method p| , the flat histogram method || 7\ , and 
the Wang-Landau algorithm H are examples of the first 
category. The second category includes the cluster algo- 
rithm; one flips a large number of spins in a correlated 
cluster at a time instead of a single-spin flip, which helps 
the relaxation time decrease drastically. Examples of the 
second category are the Swendsen-Wang (SW) algorithm 
H and the Wolff algorithm Recently Tomita and 

Okabe Jll| proposed an effective cluster algorithm, which 
is called the probability-changing cluster algorithm, of 
tuning the critical point automatically. 

The combination of approaches of two categories is 
a challenging problem to explore an efficient algorithm. 
Janke and Kappler Q proposed a trial to combine 
the multicanonical method and the cluster algorithm; 
their method is called the multibondic ensemble method. 
Quite recently, Yamaguchi and Kawashima |l3[ ] have im- 
proved the multibondic ensemble method; they have also 
shown that the combination of the Wang-Landau algo- 
rithm and the improved multibondic ensemble method 
yields much better statistics compared to the original 
multibondic ensemble method by Janke and Kappler |L2| . 

One calculates the energy density of states (DOS) 
g(E) in the multicanonical method ||, |^| and the Wang- 
Landau method ||; the energy histogram H(E) is 
checked during the Monte Carlo process. In contrast, 
the DOS for bond number rib, Ql^fe); is calculated in the 
multibondic ensemble method |12| or the improved multi- 
bondic ensemble method by Yamaguchi and Kawashima 
|jl3|| ; the histogram for bond number, H(rib), is checked 
in the Monte Carlo process. 

In proposing the broad histogram method, Oliveira et 
al. Q paid attention to the number of potential moves, or 
the number of the possible energy change, N(S, E — > E'), 



for a given state S. The total number of moves is 
N(S, E -» E + AE) = N 

AE 

for a single-spin flip process, where N is the number of 
spins. The energy DOS is related to the number of po- 
tential moves as 



E' ' 



g{E) (N{S, E - E')) E = g(E') (N(S', E' - E)) 

(1) 

where (• • ■) E denotes the microcanonical average with 
fixed E. This relation is shown to be valid on general 
grounds juj, and hereafter we call Eq. ([[]) as the broad 
histogram relation (BHR) for the energy. One may use 
the number of potential moves N(S, E — > E') for the 
probability of updating states. While the original dynam- 
ics || was criticized to be not entirely correct || |l5j , a re- 
fined dynamics is employed in the flat histogram method 
|Q . Alternatively, one may employ other dynamics which 
has no relation to N(S, E — ► E'Y but Eq. (|l|) is used 
when calculating the energy DOS |16|, [l7|]. It was stressed 
p6| , p"7| ] that N(S,E — > E') is a macroscopic quantity, 
which is the advantage of using the number of potential 
moves. We do not have to care about the relative num- 
ber of visits for different energy level E. It is contrary 
to the case of the multicanonical method EL |) or the 
Wang-Landau method ^j. The only crucial point is the 
uniformity of visits within the same energy level [p"(| . 

It is quite interesting to ask whether there is a relation 
similar to the BHR, Eq. (Q), for the bond number. In 
this paper, using the cluster (graph) representation, we 
derive the BHR for the bond number. We propose a 
dynamics based on the number of potential moves for the 
bond number. Using the DOS for the bond number thus 
obtained, we calculate the specific heat for model spin 
systems. We also employ other dynamics, that is, the 
multibondic ensemble method fl^| and its improvement 
p3[ , and calculate the bond- number DOS based on the 
BHR for the bond number. Comparing the efficiency 
of several methods, we show that the calculation of the 
bond-number DOS through the BHR gives much better 
statistics compared to the direct calculation of the DOS. 

The rest of the paper is organized as follows. In Sec. 
II, we briefly review the cluster (graph) representation for 
the Q-state Potts model. In Sec. Ill, we derive the BHR 
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for the bond number. A dynamics based on the number 
of potential moves for the bond number is discussed in 
Sec. IV. In Sec. V, calculating the accuracy of the spe- 
cific heat for the two-dimensional (2D) Ising model, we 
compare the efficiency of several methods. The summary 
and discussions are given in Sec. VI. 



II. CLUSTER FORMALISM 

We briefly review the cluster (graph) formalism for the 
Q-state Potts model. We are concerned with the Hamil- 
tonian 

H = -J y^o-i.g,) <Ji = {!,■■■ ,Q} 

(i,3) 

where J is the exchange coupling constant and the sum- 
mation is taken over the nearest-neighbor pairs 
From now on, we represent the energy in units of J, and 
the Boltzmann constant is set to be one. 

The partition function for a given temperature T is 
expressed as 

Z(T) = J2 W o(S) =J29(E)W (E(S),T) 

S E 

with the Boltzmann weight of state S having the energy 
E, 

W (S) = W (E(S),T)^e- E ^/ T , 
and the energy DOS, 

9(E) = 1 - 

{S\E(S)=E} 

With the framework of the dual algorithm [H| , the 
partition function is also expressed in the double sum- 
mation over state S and graph G as 



Z(T) = Y,Vo(G)A(S, G), 



S,G 

where A(S, G) is a function that takes the value one when 
S is compatible to G and takes the value zero otherwise. 
A graph consists of a set of bonds. The weight for graph 
G, Vo(G), is defined as 

V a (G) = V (n b (G),T) = (e 1 ^ - l) n >™ 

for the Q-state Potts model, where n b (G) is the num- 
ber of "active" bonds in G. This is nothing but the 
Fortuin-Kasteleyn representation [^o) for the Q-state 
Potts model. We say a pair (i,j) is satisfied if <7j = cr 3 , 
and unsatisfied otherwise. Satisfied pairs become active 
with a probability p = 1 — e -1 / T for given T. 

By taking the summation over S and G with fixing 
the number of bonds n bl the expression for the partition 
function becomes 

N B 

Z(T)= ]T Q(n b )V (n b ,T), 

n b =0 



where Nb is the total number of nearest-neighbor pairs 
in the whole system. Here, £l(n b ) is the DOS for the bond 
number defined as the number of consistent combinations 
of graphs and states such that the graph consists of n b 
bonds; 

n(n b )^ ]T £a(s,g). 

{G\n b (G)=n b } S 

Then, the canonical average of a quantity A is calcu- 
lated by 

(A) T - ^ , (2) 

where (A) is the microcanonical average with the fixed 
bond number n b for the quantity A defined as 

mi E {GMG)=nb} EsMS,G)A(S,G) 

Thus, if we obtain Q(n b ) and (• • -) n during the simula- 
tion process, we can calculate the canonical average of 
any quantity. 

We should note that for the calculation of the energy 
E, it is convenient to use the relation 

{E) T =T 2 ^logZ(T) = --^- J (n b ) T . (4) 
Similarly, the specific heat per one site C is given by 

,1/T 



CNT " - - (e i/r_ 1)2 <"»>t 

,1/T \ 2 



+ «"2>r -(»*&)• (5) 

The above equations (|J) and (||) were derived by Janke 
and Kappler 



III. BHR FOR THE BOND NUMBER 

The relation between the energy DOS and the number 
of potential moves for energy, the BHR for the energy, 
was rigorously derived by Oliveira Here we follow 

a method similar to that used by Oliveira to derive the 
BHR for the bond number. Instead of using the relation 
between states, we consider the relation between graphs. 

The number of potential moves from the graph 
with the bond number n b to the graph with n b + 1, 
N(S, G, n b — > n b + 1), for fixed S is equal to that of the 
number of potential moves from the graph with n b + 1 
to that with n b , N(S,G' ,n b + 1 — > n b ). That is, the 
following relation is satisfied: 

N(S,G,n b ^n b + l) = 

{G\n b (G)=n b } 

N(S,G',n b + l^n b ). 

{G'\n b (G')=n b + l} 

(6) 
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Taking a summation over states S and using the defini- 
tion of the microcanonical average with the fixed bond 
number n b , Eq. (||), we rewrite Eq. (0) as 



Q{n b ) (N(G, n b ^n b + l)) nb 
= n{n b + l) (N(G',n b + l 



n b )) nb+1 ■ (7) 



This is the BHR for the bond number. It should be noted 
that N(G, n b — > n b + 1) is a possible number of bonds to 
add, and related to the number of satisfied pairs for the 
given state S, 

rip(S) = S, 



by 



(id) 



N(G, n b — > n b + 1) = n p (S) - n b . 



With use of the microcanonical average with fixed bond 



number for n p , we have the relation 



(N(G, n b ^n b + 1))^ = {n p ) nb - n b . 



(8) 



On the other hand, the possible number of bonds to 
delete, N(G',n b + 1 — > n b ), is simply given by n b + 1, 
that is, 

{N{G',n b + \^n b )) nb+1 =n b + l. (9) 
From the BHR for the bond number, Eq. (M), we have 



fi(0) 



Hi,- 1 

n 



Z=0 



+ 1) 



n 



(JV(G,l-f l + l)> 



|= 1 (iV(G,l + l-!))„ F(+1 

(10) 

Then, substituting Eqs. @ and (0) into Eq. (|l0|), we 
obtain the bond-number DOS, fi(n&), as 

lnO(nfe) = mfi(0) + ^ In - ). (11) 



When calculating the bond-number DOS from the BHR 
for the bond number, we only need the information on 
( n p)n > ^ e microcanonical average with fixed n b of the 
number of satisfied pairs n p . It is much simpler than the 
case of the BHR formulation for the energy DOS. 

Moreover, in the computation of n p , we can use an 
improved estimator. If a pair of sites belong to the 
different cluster, this pair is satisfied with a probability 
of 1/Q. If a pair of sites belong to the same cluster, 
this pair is always satisfied. Then, we can employ an 
improved estimator Tip clS 

^^'s)?^ 1 ? 1 (12) 

where Cj(G) represent a cluster that a site i belongs to. 
Only the information on graph is needed. By definition, 



(n p ) nb = ( n p) nb - We employ the improved estimator 
in the whole calculation below. Inserting Eq. ( |T^ ) into 
Eq. (jEy), we have 



ft(n b ) 



z=o 



(1 



;)(*c(G),c,(G)) n . = , + -rt--J 



Here we have used the relation 

fi(O) = Z(T^ oo) = Q w . 



(13) 



It is interesting to check Eq. (|R|) for a special case. The 
Q — > 1 limit of the Q-state Potts model is the bond per- 
colation problem. If we substitute Q = 1 into Eq. (jllj ) , 
we obtain 



0(n b ) 



Nb 



which is the expected relation for the bond percolation 
problem. 



IV. FLAT HISTOGRAM METHOD FOR THE 
BOND NUMBER 

Let us consider the update process for the Monte 
Carlo simulation. In the multibondic ensemble method, 
a graph is updated by adding or deleting a bond for a 
satisfied pair of sites |l2]. The histogram H(n b ) becomes 
flat if we use the following rule. If there is a bond already 
on the chosen pair, we delete it with a probability 



P(n b —> n b - 1) 



tt{n b ) 



fl(n b -i) + Q,(n b ) 



(14) 



On the other hand, if there is no bond and if the pair is 
satisfied, we add a bond with a probability 



P{n b -> n b + 1) 



fl(n b ) 



Q(n b+ i) + 0(n&) 



(15) 



Since the exact form of the bond-number DOS £l(n b ) is 
not known a priori, we renew Q(n b ) iteratively in the 



Monte Carlo process by several ways 113, 13 



We may use the number of potential move for the bond 
number, (N(G, ■ ■ •))„ , for the probability of update. In- 
serting Eqs. (§),§), and (§) into Eqs. @ and @, we 
get the probability to delete a bond, 



P(n b — > n b - 1) 



( n p)n b -i + 1 - 



and the probability to add a bond, 

716 + 1 



P(n b — > n h + 1) = 



(16) 



(17) 



respectively. 
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The actual Monte Carlo procedure is as follows. We 
start from some state (spin configuration) S, and an ar- 
bitrary graph G consistent with it. We add or delete a 
bond of satisfied pairs with the probability ( |l6| ) or (f 7). 



After making such a process as many as the number of 
total pairs, Nb, we flip every cluster with the probability 
1/2. And we repeat the process. Since we do not know 
the exact form of (np) n , we use the accumulated average 
for (rip) n . The dynamics proposed here can be regarded 
as the fiat histogram method for the bond number, which 
we call the cluster-flip flat histogram method. The con- 
ventional flat histogram method for the energy will be 
referred to as the single-spin-flip flat histogram method 
hereafter. As (n p ) n converges to the exact value, the 
histogram H(n,b) becomes flat. We calculate the bond- 
number DOS by using Eq. (|l3|), and then calculate vari- 
ous quantities by Eq. (||), or Eqs.(@) and (|j). 

Here, we have described the procedure for the multiple 
cluster update of the Swendsen-Wang type || , but we can 
also employ the single cluster update of the Wolff type 
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V. RESULTS 

First, we simulate the L x L Ising model on the square 
lattice with the periodic boundary conditions by us- 
ing the cluster-flip flat histogram method. We show 
( n p)n b l-^B as a function of rib for L — 32 by the solid 
line in Fig. 0(a); we give nb/Ns by the dotted line. The 
number of Monte Carlo sweeps (MCS) is 5 x 10 7 . The 
difference between the solid and dotted lines represents 
the number of potential moves (N(rib — > rib + 1)} /Nb, 
whereas the difference between the dotted line and the 
horizontal axis represents (N(rib — > nb — 1)) /Nb- We 
should note that {n p ) n _ /Nb = 1/2, which is expected 
from Eq. (p^). The logarithm of the bond-number DOS, 
lnf2(rib), obtained by (n p ) n is shown in Fig. 0(b) as a 
function of nb- The temperature dependence of the spe- 
cific heat calculated using Eq. (||) is shown in Fig. 0; the 
deviation from the exact result obtained by Beale |2M is 
not visible in this scale. 

Let us compare the performance of the cluster-flip flat 
histogram method proposed in this paper with that of 
the single-spin-flip flat histogram method ||. To do this, 
we check the number of MCS to satisfy the flatness con- 
dition for the histogram H(rib) or H(E); we state that 
the flatness condition is fulfilled if the histogram H(rib) 
or H(E) for all possible rif, or E is equal to or larger than 
80% of the average histogram H. In Fig. ||, we show the 
size dependence of the number of MCS to satisfy the flat- 
ness condition, which we call the flatness time ifl a t here- 
after, for both the cluster-flip flat histogram method and 
the single-spin-flip flat histogram method in logarithmic 
scale. The linear system sizes L are 4, 8, 12, 16, 20, 24, 
and 32. The average is taken over many samples. The 
number of samples ranges from 20 for the largest system 
to 1000 for the smallest. We see from Fig. |^ that for the 



FIG. 1: (a) {n p ) nb /Nb and (b) lnn(n b ) of the 32 x 32 Ising 
model obtained by the cluster-flip flat histogram method. The 
dotted line in (a) denotes nt/NB- 




FIG. 2: Specific heat per a site of the 2D Ising model for L 
= 32 obtained by the cluster-flip flat histogram method. 



single-spin-flip flat histogram method the flatness time 
increases more rapidly as the system size increases. The 
least-squares fitting of the data gives 

lntfl at ~ 4.04(2) + 1.75(1) x In N 

for the cluster flat histogram method, and 

lniflat - 1.28(7) + 2.46(1) x In TV 

for the single-spin-flip flat histogram method. 

As another example, we simulate the 2D 10-state Potts 
model on the square lattice. A strong first-order phase 
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FIG. 3: Size dependence of the flatness time for the 2D Ising 
model. The linear system sizes L are 4, 8, 12, 16, 20, 24, and 
32; N — L 2 . The cluster-flip flat histogram method and the 
single-spin-flip flat histogram method are compared. 



transition occurs in this model. We show {n p ) rib /Nb for 
the 32 x 32 lattice by the solid line in Fig. ||(a); we give 
rib by the dotted line. The number of MCS is 5 x 10 7 . 
The number of potential moves (iV(n& — ► n& + 1)) /Nb 
and (N(rib —> nb — 1)) /Nb are given in the same man- 
ner as the case of the Ising model. It is to be noted that 
{n p ) nb=Q /N B = 1/10 for the 10-state Potts model. The 
logarithm of the bond-number DOS, ln£l(nfc), obtained 
by (n p ) n is shown in Fig. Q(b). The temperature de- 
pendence of the energy obtained by Eq. (Q) is given in 
Fig. [| The latent heat AQ is shown in the figure. The 
comparison of the flatness time for the 2D 10-state Potts 
model is shown in Fig. ^. The linear system sizes L are 4, 
8, 12, 16, 20, and 24. The number of samples to take the 
average ranges from 5 for the largest system to 1000 for 
the smallest. The flatness time of the single-spin-flip flat 
histogram method increases more rapidly with size than 
that of the cluster-flip flat histogram method, although 
it is not clear whether the size dependence is linear or 
not in logarithmic scale. It again shows the superiority 
of the cluster-flip flat histogram method over the single- 
spin-flip flat histogram method. 

In the calculations presented above, we have used the 
number of potential moves both for the dynamics and the 
estimator of Q(rib) or (n p ) nb . However, our procedure to 
calculate the bond-number DOS £l(rib) using the num- 
ber of potential moves, or more explicitly, using (n p ) nb , 
Eq. (jll]) or (p"3|), is independent of the dynamics. We 
may use the multibondic ensemble method Jl2| or its im- 
provement |l3"l , and monitor (n p ) n to compute il(nb), 
although fl(nb) is directly used for the probability to up- 
date and renewed with the help of the histogram H{iib), 
such as n old (nb)H(n b ) n ncw (n b ). We compare the 
accuracy of the calculation for several dynamics and the 
procedure to calculate £l(rib). For that purpose, we study 
the errors of the specific heat for the 2D Ising model. The 
energy DOS is exactly calculated by Beale pi]]. As al- 
ready shown in Fig. 0, the errors of our calculation are 




FIG. 4: (a) {n v ) nb /N B and (b) lnQ(n b ) of the 32 x 32 10- 
state Potts model obtained by the cluster-flip flat histogram 
method. The dotted line in (a) denotes n^/iVs. 







<E> T 

-1 



-2 





I 1 1 1 1 1 1 1 1 

Latent heat AQ 


J 

I 





FIG. 5: Energy of the 2D 10-state Potts model for L = 32 
obtained by the cluster-flip flat histogram method. 



very small; we treat the relative error, which is defined 

as 



e(T) 



^simulation (T^ Coxact(^') 



Cexact (T) 



for the specific heat C. The relative errors e(T) of the 
32 x 32 Ising model in the case of the cluster-flip fiat 
histogram method are shown in Fig. f^(a). The number 
of MCS is 20000 x N B - The average value of e(T) in the 
range of 1.0 < T < 4.0, which will be denoted by e(T), 
is as small as 0.0002. 

In the case of the multibondic ensemble method, we 
can calculate fl(nb) either through the number of poten- 
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FIG. 6: Size dependence of the flatness time for the 2D 10- 
state Potts model. The linear system sizes L are 4, 8, 12, 16, 
20, and 24; N — L 2 . The cluster-flip flat histogram method 
and the single-spin-flip flat histogram method are compared. 




FIG. 7: Relative errors of the specific heat for the 32 x 32 
Ising model; (a) the cluster-flip flat histogram method, (b) the 
multibondic ensemble method, and (c) the improved multi- 
bondic ensemble method. The number of MCS is 20000 x Nb ■ 
The solid line denotes the data obtained by the calculation 
using the number of potential moves, and dotted line denotes 
those obtained by the direct calculation with H(rib)- 



TABLE I: Average relative error of the specific heat e(T) for 
the 2D 32 x 32 Ising model. We compare the data for several 
Monte Carlo methods and the procedure to calculate S}(nt), 
the calculation using the number of potential moves (potential 
move) and the direct calculation with H(nt) (direct). 



e(T) 


potential move 


direct 


cluster-flip flat histogram 


0.0002 




multibondic 


0.0002 


0.043 


improve multibondic 


0.0002 


0.0087 



tial moves or by the direct calculation with the help of 
the histogram H(rib). The errors e(T) of the 32 x 32 
Ising model in the case of the multibondic ensemble 
method are plotted in Fig. 0(b). The number of MCS 
is 20000 x Nb; we renew fi(n&) for the probability of 
graph update by every 100 x N B MCS. The solid line 
denotes the data for the calculation using the number of 
potential moves, and dotted lines denotes those for the 
direct calculation using H(rib). We see that the calcu- 
lation of £l(nb) through the number of potential moves 
gives much smaller errors. The average value e(T) is 
0.0002 for the calculation using the number of potential 
moves, whereas that for the direct calculation with H(rib) 
is 0.043. We also show the results of the improved multi- 
bondic method in Fig. f?](c). The conditions are the same 
as th ose for the multibondic method. The average value 
e(T) for the calculation using the number of potential 
moves is 0.0002, whereas that for the direct calculation 
with H(ni>) is 0.0087. The calculation of f2(fi£,) through 
the number of potential moves again gives much smaller 
errors compared to the direct calculation with H(rib). It 
is interesting to notice that e(T) take almost the same 
value for several methods if we follow the procedure to 
calculate f2(nf,) through {n p ) nb - The data of e(T) for sev- 
eral methods are tabulated in Table □ for convenience. 



VI. SUMMARY AND DISCUSSIONS 

To summarize, we have derived the rigorous BHR for 
the bond number, investigating the cluster (graph) rep- 
resentation of the spin models. We have shown that 
the bond-number DOS Q(rib) can be calculated in terms 
of (n p ) nj . We have proposed a Monte Carlo dynamics 
based on the number of potential moves for the bond 
number, which is regarded as the flat histogram method 
for the bond number. We have shown the efficiency of 
the BHR for the bond number in calculating the bond- 
number DOS and other physical quantities. 

For the dynamics, the combination of the Wang- 
Landau idea ]8| and the cluster algorithms is useful in 
accelerating the diffusion of the random walker, as was 
pointed out before However, here we have made 

more emphasis on the use of the BHR for the estimator 
of fl(rib)- The advantage of using the BHR may be at- 
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tributed to the fact that the number of potential moves 
is a macroscopic quantity, which is the same situation as 
the BHR for the energy [l7| . Moreover, the use of the 
improved estimator for calculating the number of poten- 
tial moves, Eq. (|l3|), gives much better statistics for the 
calculation. 

The number of potential moves for the energy, 
N(S, E -> E±AE), has several possibilities for AE. On 
the contrary, in the case of the number of potential moves 
for the bond number, N(G, — * rib i 1), the change of 
the bond number is limited to one, which makes the cal- 
culation of the bond-number DOS through the number 
of potential moves much simpler than that of the energy 
DOS. 

Recently, a cluster Monte Carlo algorithm to simulate 
the Q-state Potts model for any real Q (> 0) was pro- 
posed by Gliozzi [^2| . It is interesting to apply the BHR 
to that method. Since only the information on graph is 



used in that Monte Carlo algorithm, Eq. (|l2|) is useful 
for calculating (n p ) . 

In this paper, we argued the BHR for the bond number. 
We can extend the present idea to the relation including 
two variables, for example, the bond number and the 
cluster number. The extension to more general cases, 
such as the loop algorithm of the quantum Monte Carlo 
simulation, may attract much attention, which will be 
studied in near future. 
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